/*
 * To change this template, choose Tools | Templates
 * and open the template in the editor.
 */

package MeshlessShapeFunctions;

import HashMatrix2.HashMatrix2;
import HashMatrix2.LUFactorization;
import Integration.IFunction1D;
import Integration.derivavty;

import java.util.ArrayList;

import jmeshless.NewOkCancelDialog;

import org.jfree.data.xy.XYSeries;
import org.jfree.data.xy.XYSeriesCollection;

/**
 *
 * @author mark
 */
public class JMeshlessInterpolationFunction1D2 implements IFunction1D {

    protected int PolinomialDegree;
    protected double xii;
    protected double rho;
    protected IKernelFunction func;
    protected HashMatrix2 A;
    protected double[] b;
    protected ArrayList<Particle> NodeList;
    protected double EPS = 0.0000001;


    public JMeshlessInterpolationFunction1D2(int PolinomialDegree, double xi, double rho, IKernelFunction func) {
        this.PolinomialDegree = PolinomialDegree;
        this.xii = xi;
        this.rho = rho;
        this.func = func;
    }
    
    public double getRho() {
		return rho;
	}

    public double getXii() {
		return xii;
	}
    
    public void setNodeList(ArrayList<Particle> NodeList) {
        this.NodeList = NodeList;
    }

    public ArrayList<Double> getBaseVector(double x, double xi){
        ArrayList<Double> L = new ArrayList<Double>();
        for (int i = 0; i <= this.PolinomialDegree; i++) {
            L.add(Math.pow((x-xi)/rho, i));
        }
        return L;
    }

    public ArrayList<Double> getBaseVectorDerivaty(double x, double xi){
        ArrayList<Double> L = new ArrayList<Double>();
        L.add(0.0);
        for (int i = 1; i <= this.PolinomialDegree; i++) {
            L.add(i * Math.pow((x-xi)/rho, i-1));
        }
        return L;
    }


    public double calculate(double x){
    	if((xii-rho>x)||(xii+rho<x)){
            return  0;
        }
        ArrayList<Double> BaseVector = getBaseVector(x, xii);
        int n = BaseVector.size();
        b = new double[n];
        A = new HashMatrix2(HashMatrix2.ARRAY_TYPE, n);
        for (int i = 0; i < NodeList.size(); i++) {
            Particle particle = NodeList.get(i);
            calculateAlfa(x,particle.x);
        }        
        double val = func.evaluate((x-xii)/rho);

        boolean key=false;

        b[0] = val;
        if(Math.abs(b[0])>EPS){key = true;}

        double returnValue = 0.0;
        if(key){
            LUFactorization lu = new LUFactorization(A);
            double[] res = lu.calculate(b);

            
            BaseVector = getBaseVector(x, xii);
            for (int i = 0; i < b.length; i++) {
                returnValue+=BaseVector.get(i)*res[i];
            }
        }
        return returnValue;
    }


    private void calculateAlfa(double x, double xi){
        ArrayList<Double> BaseVector = getBaseVector(x, xi);
        int n = BaseVector.size();
        for (int i = 0; i < n; i++) {
            double val = BaseVector.get(i)*func.evaluate((x-xi)/rho);
            for (int j = 0; j < n; j++) {
                A.addValue(i, j, val*BaseVector.get(j));
            }
        }
        //for (int i = 0; i < b.length; i++) {
        //    b[i]+= BaseVector.get(i)*func.evaluate((xi-x)/rho)*x*x;
        //}
    }

    public double calculateDiffuseDerivaty(double x){
        ArrayList<Double> BaseVector = getBaseVector(x, xii);
        int n = BaseVector.size();
        b = new double[n];
        A = new HashMatrix2(HashMatrix2.ARRAY_TYPE, n);
        for (int i = 0; i < NodeList.size(); i++) {
            Particle particle = NodeList.get(i);
            calculateAlfa(x,particle.x);
        }

        boolean key=false;
        for (int i = 0; i < b.length; i++) {
            b[i] = BaseVector.get(i)*func.evaluate((xii-x)/rho);
            if(b[i]!=0.0){key = true;}
        }

        double returnValue = 0.0;
        if(key){
            LUFactorization lu = new LUFactorization(A);
            double[] res = lu.calculate(b);

            BaseVector = getBaseVectorDerivaty(x, xii);

            for (int i = 0; i < b.length; i++) {
                returnValue+=BaseVector.get(i)*res[i];
            }
        }
        return returnValue;
    }

    public double calculateDerivaty(double x){
    	if((xii-rho>x)||(xii+rho<x)){
            return  0;
        }
        ArrayList<Double> BaseVector = getBaseVector(x, xii);
        int n = BaseVector.size();
        b = new double[n];
        A = new HashMatrix2(HashMatrix2.ARRAY_TYPE, n);
        for (int i = 0; i < NodeList.size(); i++) {
            Particle particle = NodeList.get(i);
            calculateAlfa(x,particle.x);
        }
        
        double val = func.evaluate((x-xii)/rho);
        BaseVector = getBaseVector(xii, xii);
        boolean key=false;
        for (int i = 0; i < b.length; i++) {
            b[i] = BaseVector.get(i)*val;
            if(Math.abs(b[i])>EPS){key = true;}
        }

        double returnValue = 0.0;
        if(key){
            LUFactorization lu = new LUFactorization(A);
            double[] res = lu.calculate(b);
            
            val = func.evaluateDerivaty((x-xii)/rho);
            key=false;            
            double[] tmp = new double[b.length];
        	for (int j = 0; j < this.NodeList.size(); j++) { // по вузлах
        		Particle xj = NodeList.get(j);
        			val = func.evaluateDerivaty((x-xj.x)/rho);
        			BaseVector = this.getBaseVector(x, xj.x);
        			double tmpVal = 0.0;
        			for (int k = 0; k < res.length; k++) {// альфа на вектор Р і вагу ...
        				tmpVal += BaseVector.get(k)*res[k];
        			}
        			
        			for (int k = 0; k < res.length; k++) {//число на вектор Р 
        				tmp[k]+=BaseVector.get(k)*tmpVal*val;
        			}
			}
        	val = func.evaluateDerivaty((x-xii)/rho);
        	BaseVector = this.getBaseVector(xii, xii);
            for (int i = 0; i < b.length; i++) {
                b[i] = BaseVector.get(i)*val - tmp[i];
                if(Math.abs(b[i])>EPS){key = true;}
            }
            if(key){
            	double[] res2 = lu.calculate(b);
            	BaseVector = this.getBaseVectorDerivaty(x, xii);
            	for (int j = 0; j < res2.length; j++) {
					returnValue+= BaseVector.get(j)*res[j];
				}
            	BaseVector = this.getBaseVector(x, xii);
            	for (int j = 0; j < res2.length; j++) {
					returnValue+= BaseVector.get(j)*res2[j];
				}
            }
            
            
        }
        return returnValue;
    }

    public static void main(String[] args){


        double a = 0.0,b=1;
        double x = a;


        int nn = 5;
        ArrayList<Particle> ls1D = new ArrayList<Particle>();
        double h = (b-a)/(double)nn;
        for (int i = 0; i <=nn; i++) {
            ls1D.add(new Particle(a+h*i, 0, i));
        }
//        ls1D.add(new Particle(0.5, 0, nn+1));

        ArrayList<JMeshlessInterpolationFunction1D2> ls = new ArrayList<JMeshlessInterpolationFunction1D2>();
        double rho = 1.33;
        for (int i = 0; i < ls1D.size(); i++) {
            Particle particle = ls1D.get(i);
            rho = 2.1*h;
            JMeshlessInterpolationFunction1D2 f = new JMeshlessInterpolationFunction1D2(1, particle.x, rho, new JCubicSpline(rho));
//            JMeshlessInterpolationFunction1D2 f = new JMeshlessInterpolationFunction1D2(1, 0.5, (rho)*h, new JCubicSpline(rho*h));
//            rho+=0.5;
            ls.add(f);
        }

        for (int i = 0; i < ls1D.size(); i++) {
            JMeshlessInterpolationFunction1D2 f = ls.get(i);
            f.setNodeList(ls1D);
        }

        XYSeries series1 = new XYSeries("full derivaty calculation", true, true);
        XYSeries series2 = new XYSeries("just function", true, true);
        XYSeries series3 = new XYSeries("differences", true, true);
        XYSeries series4 = new XYSeries("differences", true, true);
        XYSeries series5 = new XYSeries("differences", true, true);
        XYSeries series6 = new XYSeries("differences", true, true);
        
        
        JMeshlessInterpolationFunction1D2 jMeshlessInterpolationFunction1D = ls.get(0);        
        
        x=a;
        
        while(x<b+0.001){
        	double y = jMeshlessInterpolationFunction1D.calculateDerivaty(x);
        	series1.add(x,y);
        	series4.add(x, derivavty.calculate(x, jMeshlessInterpolationFunction1D));
            x+=0.01;
        }
        
        x=a;
        jMeshlessInterpolationFunction1D = ls.get(1);
        while(x<b+0.001){
        	double y = jMeshlessInterpolationFunction1D.calculateDerivaty(x);
        	series3.add(x,y);
        	series5.add(x, derivavty.calculate(x, jMeshlessInterpolationFunction1D));
            x+=0.01;
        }
        
        
        jMeshlessInterpolationFunction1D = ls.get(2);  
        x=a;
        System.out.println("f([");
        while(x<b+0.001){
        	double y = jMeshlessInterpolationFunction1D.calculateDerivaty(x);
        	series2.add(x,y);
        	series6.add(x, derivavty.calculate(x, jMeshlessInterpolationFunction1D));
            System.out.println(x+"\t\t"+(y));
            x+=0.01;
        }
        System.out.println("])");
//        x=a;
//        System.out.println("f([");
//        while(x<b+0.001){
//        	double y = jMeshlessInterpolationFunction1D.calculate(x);
//        	series3.add(x,y);
//            System.out.println(x+"\t\t"+(y));
//            x+=0.01;
//        }
//        System.out.println("])");
//        series1 = null;
        XYSeriesCollection col  = new XYSeriesCollection(series1);
        col.addSeries(series2);
        col.addSeries(series3);
        col.addSeries(series4);
        col.addSeries(series5);
        col.addSeries(series6);
        NewOkCancelDialog dialog = new NewOkCancelDialog(null, false, col);
        dialog.setLocationRelativeTo(null);
        dialog.setVisible(true);

    }

}
